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1. Introduction 


This report Is the annual report on work done under NASA Grant NSG-1603, 
NASA Langley Research Center, Hampton, Virginia. The grant task Is to In- 
vestigate nvunerlcal algorithms for analysis, and design of large space 
structures. The report covers activities during the period March 16, 1980 
to March 1^^, 1981. The report for the previous period, March 16, 1979 to 
March 15, 1980 was on matrix polynomials and their relationship to matrices 
of the form sI-A. The report during the more recent period uses that material 
to derive new results. 

Section 2 of this report describes the sign algorithm and the application 
to decoupling of differential equations such as the finite element equation 

(1.1) tfic(t)+Cx(t)+Kx(t) = u(t) 


The algorithm Is useful for the decoupling procedure but does not appear to 
be efficient for large scale structures. Work Is continuing on finding new 
approaches to the decoupling (or model reduction) problem. 

The work given In Section 2 Is extended In Section 3 with new results. 
The generalized sign algorithm Is given and Its application to several prob- 
lems discussed. 

Sections 4 and 5 discusses the Laplace transforms of matrix functions 
and the dlagonallzatlon procedure for the equation given In (1.1). It Is 
shown that the matrix 


( 1 . 2 ) 



I 


-C 


J 




1 


2 


can be transformed to the form 


(1.3) 


BD 


0 I 

- 2^03 


by a constant matrix Q. The matrix polynomial requires a matrix polynomial 
to diagonalize with 


(1.4) 


(Is^+2^o)s-Ho^) 


Q(s) [ls^+Cs+K]P(s) 


Extensions of the results In Section 5 will be studied during the period of 
the present grant. 

Identification of the mass, damping and stiffness matrices, M, C and K, 
la considered In Sections 6 and 7. It Is shown that the quadrature Inte- 
gration algorithm can be used to determlt.e the elements of the matrices. 

It would appear that the algorithm Is suitable for large space structures 
although the algorithm requires further examination. Actual data will be 
used with the algorithm to determine the accuracy of the procedure. 

Some results have been obtained on the computation of the damping matrix 
C for prescribed eigenvalue locations. The work on this matter Is not re- 
ported due to the timing of the report as well as the need for additional 
research. This problem V7lll also be explored during the remaining period 
of the present grant. 

An Interim report was filed earlier on the Ibrahamln Identification 
algorithm. Readers Interested In that work should secure a copy from NASA 
Langley Research Center. 


2. The Sign Matrix and Decoupling 


There has been considerable Interest In the analysis of dynaulc systems 
by decoupling or order reduction, [lj-[Aj. The decoupling of a system Into 
subsystems with different time scales has been described earlier by Yoo, 

[5] and by Popeeva and Lupas, [6] where a block dlagonallzation scheme vas 
used to separate the modes . The mathematical theory of the block dlagonal- 
Izatlon procedure was given by Russell, [7] In a earlier paper. The concept 
of separating "slow" and "fast" modes was the basic motivation for the work 
by Lee, [3J, and Yackel and Kokotovlc [8] although their approach differs 
from the procedure given In [l]-[6] and the algorithm to be presented. 

This section will describe a general procedure for decoupling of a 
system Into subsystems. The number of modes In each subsystem Is arbitrary 
with the eigenvalues arranged according to their magnitudes. The resulting 
subsystems will not be stiff even though the system may be stiff. The 
system dynamics Is recovered from the subsystem d3mamlcs by using the de- 
coupling matrices. The method Is closely related to the work of Kron, [9], 
and system decomposition, [10]. The concept as given here was first given 
by Yoo, [5]. 

The algorithm to be presented la based on the solution of algebraic 
Rlccatl equations and the sign function, [4], conceived by Roberts. The 
sign function and Its application to system analysis ha.3 been Investigated 
by Beavers and Denman [11]-[15], as well as by others, [16]-[21]. It can 
be shown that the spectral domain p(A) of a system matrix can be decomposed 
Into spectral subdomains, p^(A), by the sign matrix of A with the union of 
Pj^(A) covering the domain p(A). The eigenvalues of A remain Invariant under 
the decomposition operation and the elgenprojectors [22], for the subdomains 
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can also be generated as a byproduct of the decomposition. 

The sign function for a system matrix will be described In this papers 
as well as the decomposition procedure. The concept of elgenprojectors will 
be Introduced as well as the numerical use of elgenprojectors. 

The system equation to be considered will be of the form 


( 2 . 1 ) 


M 


dt 


dx(t) 

dt 


+ K x(t) 


f(t) 


where M, C and K are the mass, damping and stiffness matrices respective- 
ly. The vectors x(t) and f(t) are nxl with M, C , and K n^n matrices. 

The system matrix will be written as 


( 2 . 2 ) 


A = 


0 I 


L-K -CJ 

where A Is 2nx2n with K ■ -M , and C - M . The system matrix A will 
have 2n eigenvalues and It will be assumed that all eigenvalues In the 
spectrum p(A) occur In conjugate pairs with an even numb. r of zero values. 

The form of system equation assumed occurs In structural dynamics. 

The procedure will be valid for A In general ^.orm, the above form was 
chosen since the procedure is applied to structural dynamics In the current 
research effort. 

Section 2.1 will describe the sign function which has an Important 
role In the decoupling algorithm. The elgenprojectors, or projectors, will 
be Introduced In Section 2.2 along with the application of elgenprojectors 
to numerical analysis. Decomposition of the spectrum of A Into subdomains 
will be discussed In Section 2.3. The application of the decoupling al- 
gorithm to system equations Is given In Section 2.4. 


2.1 The Matrix Sign Function and Projectors 


Let A be a 2nx2n matrix with 2n complex eigenvalues where Is the 
1th eigenvalue. Assume that J denotes the Jordan form of the eigenvalue 
matrix with Jordan blocks J^ for 1 ~ 1,2, . . . ,s^2n. An eigenvalue may be 
repeated with multiplicity and the A matrix may have a null space less 
than 2n requiring generalized eigenvectors where denotes the number 
of generalized eigenvectors for X^ and s Indicates the number of Jordan 
blocks of A. 

The definition of a function of a matrix Is given by Gantmacher, [21]. 
If A and J are similar matrices such that <t> transforms J Into A, 

(2.1.1) A - •I'J'I’"^ . 

Then the matrix functions f(A) and f(J) are similar and 9 transforms f(J) 
into f (A) ; 

(2.1.2) f(A) - $f(J)$"^ . 

Defining the sign of a complex function as the sign of the real part by the 
complex function, then sign (A) Is given by 

(2.1.3) S ■ slgn(A) ” $slgn(J)$ ^ 

where Is the matrix of eigenvectors of A. Since J Is the Jordan fom, 
the eigenvalues of A will be located on the diagonal of J with pl<is ones on 
the off-diagonal of J provided that such Is required for multiple elgen- 
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values. The sign of J will be defined as 

(2.1.4) sign (J) ■ sign (A) 
where A ■■ dlag[Jj thus S Is defined as 

(2.1.5) S ■ slgn(A) ■ slgn[Re(A) ] ^ 

Roberts, [4], has shown that S can be computed by the Iterative algorithm 

(2.1.6) - J IS^+(S^)"^] S° - A 

2 

which Is the Newton algorithm: for S » 1. This algorithm will have quadratic 

2 

convergence and the Iterative process can be terminated when trace (S )-n<e. 

It should be noted that the sign of Re(l^) does not exist when ■ 0 or 
X^ « Jb)^. A method of computing the sign under such an eigenvalue assignment 
will be described later. 

The elgenprojectors of a matrix are defined as the matrices 

(2.1.7) ^ “ l,2,...,s<2n 

for the S Jordan blocks and secondary elgenprojectors with 

(2.1.8) P^j - j - 1,2 

for the repeated eigenvalues X^ where Is the number of generalized eigen- 
vectors for that eigenvalue. The matrix Is given by 
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(2.1.9) - diag[0 0 ... C 1 ... 1 0 ... 0] 


uhere the ones occupy the same positions as does In the 1th Jordan block. 
The matrices are defined as the set of matrices with ones on the off- 
dlagonals of the matrix. To Illustrate, let be defined as 


( 2 . 1 . 10 ) 



then 



» M 

10 0 


0 1 o' 


0 0 1 

(2.1.11) E^q - 

0 10 

"ll- 

0 0 1 

1 

CM 

0 0 0 


0 0 1. 


0 0 0 
^ J 


o 

o 

L_ 


The E^q matrix can be expressed In the compact notation 


( 2 . 1 . 12 ) 


“ dlaglOy 0, 0 ... 0, 


with a similar expression for E^^ but where E^^ will have representations, 

^iO’^il 

The primary eigenprojectors, will have rank r^ where r^ Is the 
multiplicity of the eigenvalue. The secondary eigenprojectors, 
j ■ l,2,...,il^, will have rank ^he eigenprojectors are idempotent 




and have properties; 
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S 


2P1. 

'lO • ' 

1-1 


2P2. 

^10^10 “ **10 

(Idempocent) 

2P3. 

^io^jo ’ ° 

1 ^ J 

2P4. 

P p « p 

ij Ij 1,J+1 

j - l,2 ,...i^-l 

2 P 5 . 

P P - 0 

*^ 1 J*^ 1 ,J +1 

J “ Xi 2 ,..., 


It is not difficult to show that any square matrix A has a spectral decom- 
position 


(2.1.13) 


1“1 


A * 

*i 


s 

I 

i-1 


10 


- y 

i-1 


^ 1^10 


s 

+ I 

1-1 


il 


where is the set of primary elgenprojectors and is the first second- 
ary elgenprojector of a Jordan block if it exist. 

The elgenprojectors of A, which is assumed to be time- Invariant, can 
be computed by repeated computations of the sign matrix. Assume that m 
eigenvalues of A have eigenvalues with Re(X^)>0 and 2n-m eigenvalues have 
Re(X^)<0. This spectral splitting of the spectral domain can be achieved 
in several ways, origin shifts, bilinear transformations, etc. The sign of 
A will then have the form 


(2.1.14) 


S - $ 


mxm 
0 -I 


.-1 


2n-mx2n-m 


The lower Identity matrix will vanish if 0.5(I*fS) is computed whereas the 
upper identity block will not be present in 0.5(I-S). The projector for 




9 


Che positive and negative spectra are then given by 

(2.1.15) - I (I+S) 

(2.1.16) i • 

^ "t* * 

The products AP and AP gives A and A respectively where the nonzero 
eigenvalues of A^ have Re(X^)>0 and those of A have Re(A^)<0. The spectra 
of A and A can now be altered by an origin sh^ft and the signs computed. 
Repeated con 9 >utation of the sign matrix, origin shift, and spectral decom- 
position will give the eigenprojectors of all Jordan blocks if desired. 


2.2 Decomposition of Che Spectrum of A 


The spectrum of e general nxn matrix A can be decompoaed into m matricea 

with subspectra p . (A) such that p(A) “ P, (A) U ^2(A) . . . U P (A) where m<s. 

1 1 o — 

Each matrix A^ will have an associated set of r finite eigenvalues depending 
on the spectral decomposition as well as n-r zero eigenvalues. The decom- 
position of A is given by 


( 2 . 2 . 1 ) 


s 8 m 

A ■ I ^l^io * ^ 'll ■ ^ *1 
1-1 ^ i -1 1-1 ^ 


where can be any partial sum over a selected set of eigenvalues. The 
decomposition of A may be by eigenvalue magnitudes, reals, complex, imaginary 
or any other desired eigenvalue property. 

Suppose that all eigenvalues with magnitudes 
assigned to a spectral domain p(A^), then A can be decomposed into m matrices 
such that 


( 2 . 2 . 2 ) 


- ! A 

J-k J 


i — 1,2,... ,m 


with 


m m 

(2.2.3) A - I A - A I P. 

1-1 ^ 1-1 ^ 

Each A^ will be 2nx2n with 2n-(p-k) zero eigenvalues and p-k nonzero eigen- 
values in the spectral domain p^(A). 

The spectral decomposition described above is a decoiq> 08 ltlon of A into 
m 2nx2n matrices. It is also possible to find m square matricea of lower 
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11 


dimension that will preserve the spectral domain of A. Let the eigenvector 
matrix <I> be partitioned as 


( 2 . 2 . 4 ) 


11 

12 

L "^21 

CM 

CM 


where and ^^2 square matrices of dimension mxm and 2n-mx2n-m re- 
spectively. Define the transformation matrix T as 


( 2 . 2 . 5 ) 


I ^{diag[4’^ - *22 ]} 


1 

2 


4 $ ^ 

21 11 


*12 22 


-I 


The similarity transformations T ^AT and TAT ^ will produce block diagonal- 
ized matriceu. The first transformation T ^AT gives 


(2.2.6) Ag - t’^AT - diag[4>^j^-4>22]'J diagl*"^-*^^] 

- diag[Ag^ Ag2J 

The Ag^ submatrices are equal to where Jg^ contains the Jo<.dan 

blocks associated with the eigenvalues belonging to the set of eigenvectors 
in If has eigenvalues *^32 eigenvalues with 

|X^|<Pj^ then Agj^ has only "fast" modes with the "slow" modes In Ag 2 * 

The transformation matrix T can be constructed from the eigenprojectors 
of A. Let p'*’ denote the sum of the eigenprojectors P^^q for eigenvalues 
|Xij>Pi and P the sum of the remaining eigenprojectors. The block eigen- 
projector P"^ is given by 



12 


(2.2.7) - J PjO “ 


k 


where I - dlagll.O] * dlag[E ^,0] - dlag[ > E.-,0]. The Indexing over the 
1 Bi 10 

susanatlon assumes Chat the eigenvalues of A have been arranged by descending 
magnitudes. Since the sum of the elgenprojectors must be equal to the Iden- 
tity matrix, then 


( 2 . 2 . 8 ) 


2n 

I 

1-k+l 


PlO - - (I-P*^) 


It is not difficult to show that the sign of A can be computed from the 
elgenprojector , or the elgenprojectors can be determined from the sign ma- 
trices. Assume that k eigenvalues have Re(A^)>0 and 2m-k eigenvalues have 
Re(X^)<0. If is Che sum of the elgenprojectors for Re(X^)>0 and P the 
sum of Che elgenprojectors for Re(X^)<0 then 


(2.2.9) 


S - sign (A) - P'*’-?” 


k 2n 

io ‘■‘o ■ iL« 


where P^^ is the individual elgenprojector for a Jordan block. Noting that 
P ■ I“P^» then 


( 2 . 2 . 10 ) 


S =■ slgn(A) - 2P -I 


or 


(2.2.11) P'*' - y (S+I) 

with P~ determined from p"*". The elgenprojectors P^^ can be computed by 
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utilizing origin shifts or bilinear transformations. These procedures will 
be discussed later. 

The transformation matrix T was defined in (2.2.5) and this matrix can 
be computed from the sign matrix o'* the eigenprojectors. According to 
(2.2.9), the sign of A Is given by P^-P thus 

(2.2.12) S » P'^-P" - 4» slgnlRe(A)]4>"^ 

provided that Re(A^) ^ 0 for all 1. Let sign [Re(A)] ■ E ■ dlag[l^,-l 2 ] 
then 

(2.2.13) T = [S+E]“^ - [«1)E'I>""+Ej"^ “ ♦IJE+E'I']"^ 

since [$E+E4*] ^ - dlag[4>^^-4>22] • should be pointed out that the signs of 
(2.2.5) and (2.2.6) can be reversed. 

The analysis thus far has considered only block dlagonallzatlon for 
two subblocks. It Is not difficult to show that the procedure Is valid for 
s blocks, the number of Jordan blocks In A. The limiting case Is equivalent 

- dlagI\j^,X2 

the eigenvalues are distinct. 

Knowledge of any one of the matrices, T, S, P^ or P , Is sufficient to 
coiq>letely determine the other three matrices for two block, dlagonallzatlon. 
If m blocks are to be found, then additional Informction must be available 
from the eigenprojectors or the sign matrix. Details on the spectral de- 
composition will be found In [23]. 


to scalar dlagonallzatlon of A In which case A^ 


2.3 Matrix Sign Functions and Block Dlagonallcatlon 


It was shown In the previous section that spectral dacoivoslclon and 
block dlagonallzatlon Is possible If the elgenprojectors are known. It will 
be shown In this section how the matrix sign function Is utilised for carry- 
ing out the block dlagonallzatlon. 

Assuae that A has real and complex eigenvalues and A la to be block 
diagonalized. It will be assumed that the spectral domains are such that 
|X^|<a will be In p(Ag^) and |X^|>a will belong to p(A^ 2 ^* '*** 

strlction on the procedure Is that a be selected such that there la no 
eigenvalue with magnitude equal to a. The bilinear transformation will map 

(2.3.1) Ap - (A-al) (A+al)"^ 

all eigenvalues |X^|<a Into the left half plane of the cosq>laz apace and all 
|\^|>a Into the right half plane. The spectral domain of ^ 

entire left half plane and that of P 2 (Aj^) 3Ls the right half plane with 
p(Ajj) “ j^(Ajj) U p 2 (Ajj) . The sign of Aj^ Is then computed by (2.6) with 

(2.3.2) S » sign (Ap) - ${diag[Ij^ l 2 J)*"^ ‘ 

/N /V 

The matrices and I 2 will not be ordered In the general case with ±1 
intermixed along the diagonal of the two matrices. The reason for this la 
that the eigenvalues of A are not ordered and the Iterative algorithm for 
the sign of A^^ Is Independent of the order. 

Assuming that I^ has only plus ones which Indicates an eigenvalue |X^|>a 
and I 2 has only minus ones for |X^|<a then S can be written aa 


14 


15 


(2.3.3) S - sign (A^j) - 4>E<1> 


where E Is an elementary matrix. It follows that 


(2.3.4) S+E - <I>E$”^+E - I<I'E+E4>]4>' 



and the Inverse of S+E Is 


(2.3.5) 

T =■ [S+E]"^ - y $ 

11 

0 

^-1 

1 

2 

I -Rj^2 



0 

22 -J 


. ^21 - 


where “ ***12*22 ^12 “ *21*11 * 

The similarity transformation T ^AT will produce a block diagonalized 

matrix similar to A, [24], with 


(2.3.6) 



^ ll '*’^ 12*^21 ^ 12 '*'^ 11 ^ 12 "*^ 12 ^ 22 "^ 12 ^ 21*^12 

.^ 2 1'^^2 2 *^ 2 1 "^ 2 1 ^ 1 1 "*^ 2 1 ^ 12^2 1 ^2 1 ^ 12 _ 


where ^21 solutions to the two algebraic Rlccatl equations, 

[13], 


(2.3.7) 

(2.3.8) 


^ 12 ‘^^ 11 ^ 12 "*^ 12 ^ 22 "^ 12 ^ 21*^12 ” ° 
^2l'^^22^2l“*^21^1l"*^21^12*^21 " ° 
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The solutions to these equations are given by the off diagonal blocks of T 
as seen In (2.3.5). It Is a simple matter to show that the mode decoupling 
algorithm of Anderson, [1] Is a special case of the above analysis requiring 
a two step process, computation of R ^2 ^21^ solution to an 

associated Lypanov equation. 

The dlagonallzatlon procedure given above was for dlagonallzatlon of 
A Into two blocks with the eigenvalues | j >a In A^^ and all others In A^ 2 • 
Suppose now that the spectral domain Is to be factored Into three domains 
with ^2^^^ 1^1 1 ^“2 P 3 (A). The pro- 

cedure Is the same as described for the two block decomposition except that 
A^. Is now used In the bilinear transformation with 

(2.3.9) Ajj^ = (Ag^-a2l) (Ag^+<X2l)"^ . 

The sign of A^^ Is conq>uted and for the sign (A^j^) Is then used to find 
T^ with 

(2.3.10) - (Sj^+Ej^) 

The transformation matrix T^ Is augmented with an Identity matrix I having 
the dimensions of A,^ with 

(2.3.11) T^ - 

If Tq denotes the first transformation matrix and T^ the second, then the 
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product must be 


(2.3.12) 


T T “ — 
^0 1 4 


*"21 


-R 

-I 

-R 


12 "13 


32 


23 

I 


from which It follows that 


(2.3.13) Ag 


^ll‘^^12^2l'’'^13^31 

0 

0 


^ 22 '^^ 21 ^ 12 ''’^ 23*^32 


0 

0 


^ 33 '^^ 31 *^ 13 '^^ 32*^23 


The R^j matrices satisfy coupled algebraic Rlccatl equations, [16]. 

An ex.unple Is given to Illustrate the above procedure. Let A be de- 
fined as [25] 


A - 



-2 0-2 5 

-3 2 0 6 

8-5-6 0 

-5 5 1-2 

-6 1 6-3 

0-2-3 8 


with eigenvalues “ “1*598734, 13*1^ “ 4.455989 and “ 16.142744. 
Values of ■ 2 and 8 were selected to split the spectrtim. The block 
matrices are 


^1“ 


16.1427 

0 


0 

16.1427 
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^2 " 


4.45599 0 

0 4.45599 


^3 “ 


-1.59873 0 

0 -1.59873 


The transformation matrix for the dlagonallzatlon Is 


T 


0.25 

0 . 

1.51992 

-0.116553 

-0.678077 

-0.0496661 


0 . 

0.25 

-1.2999 

0.0515402 

0.641664 

0.379825 


0.595655 

-0.4012 

-0.25 

0 

-0.0845469 

-0.0885163 


0.407491 

-0.431213 

0 . 

-0.25 

0.214427 

0.0860365 


0.530013 

-0.407491 

-0.397965 

0.231053 

0.25 

0 . 


0.120331 

-0.37346 

-0.167369 

0.0710525 

0 . 

0.25 


which Is since A Is diagonalized. The computations were made In double 
precision. 

An Important aspect of the Integration procedure Is controllability of 
a mode. It can be shown that the system modes are controllable If and only 
If the vectors In B spans the space of A. Let B be 2nXm with vectors b^ 
such that 

(2.3.14) B “ [b. b_ ... b ] 

X z m 

The mode exp(X^t) Is controllable If the sum of the column vectors In B Is 
linear dependent on the eigenvector for To show this, consider the 

Laplace transform of z(t) ■ Az(t)-fBf(t) which is 
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8 IC c K 

(2.3.15) Z(s) " IsI-A]"Hf(s) - { I -r^ + I 

1-1 ® ®1 J-1 

where denotes the primary matrix residues 


(2.3.16) " (s-s^)’ [sI-A]"^ 


1 — lf2y...fS 


8-S , 


and the secondary matrix residues 


(2.3.17) 


■‘ij ■ jr^ 


J — If 2( . . . »r—l 


s-s. 


where r Is the multiplicity of an eigenvalue. It then follows that If 
and are zero for all j, the 1th mode will not appear In z(t). Since 

the columns of B can always be constructed from a linear combination of the 
eigenvectors of A, which must span a 2n dimensional space, and since the 
residues are elgenprojectors , the absence of the 1th eigenvector In B will 
make K^qB and K^jB zero. Modes that are not controllable will require the 
Inclusion of f(t). 
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If f(c) !• rapidly varying function with respect to a node, the integration 
step will be set accordingly. 

Since A has the special form given in ( 2 . 2 ), the first n eleaents of 
z(t) will be the displacements x(t) and the last n elements of the vector 
will be x(t). If only x(t) is desired, it Is easily shown that 

( 2 . 4 . 4 ) x(t) “2^ [v^(t)-R^2'^2^*=^‘^*^13''3^^^--'^ 

where m is the number of diagonal blocks in A^. 
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If f(c) Is rspldly varying function with respsct to a mods, tha Intagration 
step will be set accordingly. 

Since A has the special form given in (2.2), the first n elements of 
z(t) will be the displacements x(t) and the last n elements of the vector 
will be x(t). If only x(t) is desired, it is easily shown that 

(2.4.4) x(t) - 2 ^ 

where m is the number of diagonal blocks in A^. 




2.5 Application to "raa- Fra* Saao Modal 


Tha algorithm waa appllad to a frae-fraa baaa modal ganaratad by NASA 
LRC. Tha baam waa 144 Inchaa in length, tha maaa per unit length waa 0.112 

5 

Iba with tha El conatant equal to 6.25^10 . The nodes ware aqul-spacad 
with spacing with tha two end nodes located at a distance of 1/2 from 
the ends. A total of 21 nodes was assumed although the computer model 
allowed the number of nodes to be changed. 

The decoupling algorithm was used to separate the highest frequency 
mode from all others. In addition to obtaining the solution by decoupling, 
the initial conditions were chosen such that only the highest mode was ex- 
cited. The full set of equations were then integrated to obtain a benchmark 
solution for comparison to the decoupled solution. A step size of h*0.0002 
waa selected for the two computed solutions. The displacement for the hip- 
est mode is given below at three selected times where v(t) is the transformed 
solution vector v(t) - T ^z(t), as given in Section 2.4 where z(t) is the 
displacement of the free-free beam nodes. 

t V 2 ]^(t) (Full model) V 2 ^(t) (Decoupled) 

0.002 0. 584225 754E-3 0.58425754E-3 

0.004 -0.11151369E-2 -0.11151369E-2 

0.006 -0.13946460E-2 -0.13946460E-2 

The difference between the two solutions was in the 11th digit thus it can 
be assumed that the decoupling algorithm performs in a satisfactory manner. 

Test of the algorithm has been carried out for several other models and 
the accuracy was excellent. There is no reason to believe that the algorltlun 
falls for any model, the major problem with the algorithm is the computational 
overhead in computing the decoupling matrix T. Work continues on the al- 
gorithm but it appears that the algorithm is not suitable for large space 
structures in its present form. 
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(3.6.5) x(0) ■ (P^+P +P^+P^)x(0) ■ x^(0)+x (0)+x^(0)+x^(0) 


The solucion Co Che syscem equaclon is Chen given by 


(3.6.6) x(c) - exp(A C)x"*^(0)+exp(A C)x (0)+exp(A c)x^(0)+x® (0) 


where'che propercles of Che elgenprojeccors has been uClllzed Co obcain (3.6.6) 
from (3«6.4) and 0*6.5) . Equacions 0 *6.4)-0 .6.6) are Che machemaclcal expressions 
for Che decoupling concepC. Ic is possible Co conCinue Che decoupling process 
furcher such chac Che decoupled solucion vecCors include "fasc" and "slow" 

'*<odes alchough chac procedure will noc be given since iC has been described 
elsewhere, [6] and SecCion 2 of Chls reporc. 

The following maCrix has been used for llluscraclon of Che algorlchm 
where A has Che following eigenvalues ■ 1.0 ■"3.0 A^ ■ +2J A^ ■ -2j 


(3.6.7) 


0 1 0 0 ~ 

0 0 10 
0 0 0 1 

12 -8 -1 -2 


and Che IniCial condlClons for Che sysCem X ■ AX am given by X(0) ■ [-1.0, 
T 

4.0, 1.0, -3.0] . The resulCanC eigenprojecCors are 


0.6 

0.2 

0.15 

0.05 

0.6 

0.2 

0.15 

0.05 

0.6 

0.2 

0.15 

0.05 

0.6 

0.2 

0.15 

0.05 





3. G«narallz«d Sign Matrix 

Applications of tha aatrlx algn function hava baan ).'ravlously glvan In 
[15]. Tha matrix algn function la a uaaful mathematical tool for apactral 
dacooiposltl I as it provides Insight on tha spectrum of the alatrlx even 

/ 

though tha eigenvectors and eigenvalues are unknown. J 

I 

The material given In this paper Is an extension of earlier work with 
the concept of a generalized sign matrix Introduced and discussion of several 
applications. Knowledge of the generalized sign matrix Is sufficient to 
V natruct the elgenproj actors on the spectral domain of a nxn matrix. The 
generalized sign matrix appears to have been Introduced by Mattheys* [19] 
as the extended matrix sign function. The computation of the generalised 
sign mstrlx Is given in the next section. 


3.1 The Generalized Sign Function of a Matrix 
A general nxn A matrix can be represented by the form 

(3.1.1) A - MJM"^ 

where M is the eigenvector matrix of A and J is the Jordan canonical form 
with the eigenvalues along the diagonal and ones or zeros on the auperdlagonal. 
If J is expressed as 

(3.1.2) J - diag[J(Xj^)^J(X2)y.jJ(^r^l 

Then any function f(A) which is analytical in the plane containing the eigen- 
values of A is defined [3] by 

(3.1.3) f(A) - Mlf(J))M"^ 


where 


(3.1.4) f(J) = dlag[f(J^),f(J2)..,f(J^)] 


wJ’.li each Jordan block having the form 


(3.1.5) 


f(A^) f’(A^)/ll 

0 f(X^) 

• • 

• • 

• • 

0 0 


(n ,-i) 

f ^ (X^)/(n^-l)I 

(n.-2) 

f ^ (Xp/(n^-2)! 


f(X,) 
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I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


te-'v.'-B-’ir*' 


26 

Lf the sign function for a Jordan block is defined as the function with the 
sign of the real part of It then follows that 

(3.1.6) sign (A) - II sign(R^J]M ^ 


where aign [R^J] ■ diag [signCR^Xj^) sign(R^l 2 ) • • *sign(RgX^) ] . 
with " 


(3.1.7) 


sign(R^X^) 


f +1 Rg(X^)>0 

[-1 R^Up<0 


The definition of the sign of A has a restriction in the sense that the 
sign of zero or imaginary eigenvalue is not defined. To account for this 
case define the generalized sign matrix [19] as 


C3JL.8) 


r 


1 


sign(X^) 



R ^( X^)>0 

R^(Xp<0 


with the generalized sign matrix of A denoted by S(A) with 


(3.1.9) S(A) - M sign(J)M‘ 


An cfriclciiL algurithm tu find tlie sign uf a matrix lias been given by 
Roberts [4] which is based on Newton's algorithm for = 1. The iterative 
algorithm 


S 


i+1 


1 

i 


CSi-s-/) 


A 


(3.1.10) 
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I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


has quadratic convergence with convergence In approximately k-t*3 iterations 

k -1 2 

where 2 >{max and (min A^) }. The stopping criterion iTr S -nl£e has 

been used by the authors where e is a preselected error bounds. Several 
versions of accelerated convergence algorithms have been given in the liter- 
ature [4], [20], 

It should be noted that the iterative algorithm in (3.1.10) will fall 
when A has a zero or imaginary elgenvi.lue. The singularity encountered in 

(3.1.10) can be avoided by a translation of the spectrtjun or an origin shift 
by adding pi to A with 

(3.1.11) A+pI » MAM"^+PMM~^ - M(A+pI)m”^ 

wliere the direction of the shift is set by the sign of p. There is no assur- 
ance tliat a random choice of p will not produce a new singularity thus care 
should be taken in selecting p. 

Assume that p is selected such that p<|R^A^| for all eigenvalues having 
a nonzero real part and furtlier, let A+pI and A-pI be nonsingular. Let S(A^) 
and S(A 2 ) denote the signs of A+pl and A-pI respectively. The generalized 
sign then given by 

(S(A,)+S(A,)) 

(3.1.12) S(A) = ^ ^ 

The cumputationai strategy is to move all cigenv.-lucs along the Jw axis into 
tlie riglit half pleinc for S(Aj^) and into the left hand plane for S(A 2 ). The 
eigenvalues with K^(^j_) 0 must remain in the same spectral domains under 

the sliifts. The major task associated with computing the generalized sign is 
the selection of p. 

* 

I 




3.2 Eigenprojectors of the Sign Spectrum Decomposition 
The ability to compute the elgenprojector from the generalized sign matrix 
makes the sign matrix a useful tool for system analysis. The availability of 
the eigenprojectors makes It possible to generate new system matrices based 
on various mathematical operations on the eigenvalues of the system. The 
eigenvectors will remain Invariant under the operations « [22j . 

Suppose that the eigenprojectors are defined to cover the entire spectrum 
and let the eigenprojectors be denoted by P^, P , P^ and P^ with 


3.2D.1 
3. 2D. 2 p" 
3. 2D. 3 p® 


Positive elgenprojector > the projection on the positive real 
part of the eigenvalues of the matrix; R^(X^)>0. 

Negative elgenprojector - the projection on the negative real 
part of the eigenvalues of the matrix; R^(X^)<0. 

Null elgenprojector - the projection on the zero eigenvalues 
of the matrix. 


3. 2D. 4 p^ 


Imaginary elgenprojector - the projection on the Imaginary 
eigenvalues of the matrix. 


when these eigenprojectors are applied to the A matrix, a set of subspaces 
will be obtained with 


A - AP = P A 


A = AP = P A 
A® = AP° “ P°A 


A*^ = AP^ = P^A 


with A *= A^+A~+A^+A^. It should be obvious that A® is a zero matrix but 
will not be a zero matrix. The eigenvalues of the eigenprojectors will be 
zeros and ones with the same eigenvectors as A. The definition of these 
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projectors as eigenprojectors follows from the eigenvector property. The 
elgenprojectors have the basic properties 

3.2P.1 the sum of the eigenprojectors is the Identity matrix, P^+P^+P +P'*’-1 

3.2P.2 the eigenprojectors are orthogonal, i.e. P P^ •> 0 

■i^P»3 the eigenprojectors are idempotent matrices, i.e. P"*^P^ ■ P"*^. 

These 'three properties are useful in system analysis and provide the essential 
information for spectral decomposition of A. 


I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


3.3 Computation of the Elgenprojectors 
The elgenprojectors defined in Sections. 2can be computed in a straight- 
forward manner from the generalized sign matrix. The generalized sign matrix 
must be such that 

(3.3.1) S(A) - 

Since all eigenvalues along the jw axis will have sign(R^X^) •• 0. The square 
of the generalized sign matrix must be 

(3.3.2) S^(A) - P'^+P" 

Since the elgenprojector are orthogonal and Idempotent matrices. It then 
follows that P"^ and P are given by 

(3.3.3) 

(3.3.4) p- . i-t-AjJlA). 

Recognizing Chat Che sum of the elgenprojectors must give the identity matrix, 

It follows that 

(3.3.5) P^+P^ = l-p‘*’-p" = I-S^(A) 

To compute P^ and P^, the eigenvalues in the domains of P^ and P are eliminated 
by multiplying (3.3.5) by A thus 
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(3.3.6) 


^ 0 I 0 1 ^2 

- A>A - ACP'^+P ) - A(I-S^(A)) 


0 

where A ■ 0. The square of A^ will have only real eigenvalues Cherefori. 

/s 

the matrix A^ defined by 

(3.3.7) A 2 - A+aJ - A+A^(I-S^(A))^ - A+A^(I-S^(A)) 

will have a spectrum with only real, complex and zeroelgenvalues . All Imagln- 

2 

ary eigenvalues will have been shifted by The original 

eigenvalues with ^ ^ will be preserved In the operation. The eigen- 

projector for the imaginary spectral domain will then be given by 


(3.3.8) . S(A)-S(A 2 ) 


and finally 


(3.3.9) P° » I-P'''-P"-P^ - I-S(A)-S^(A)+S(A 2 ) 


The latter elgenprojector , P^, can be utilized for Inverting singular matrices 
as will be shown later. 

Example 1 - The following matrix is used to find its spectral decomposition 


which has elj'cnvalues X^ “ 0 

II 

CN 

- 2.0 X^ 

= 3.0 

^4 “ +J 


29.2 

-24.2 

69.5 

49.8 

7.0 


- 9.2 

5.2 

-18.0 

-16.8 

- 2.0 

A « 

- 10.0 

6.0 

- 20.0 

-18.0 

- 2.0 


- 9.6 

9.6 

-25.5 

-15.4 

- 2.0 


9.8 

- 4.8 

18.0 

18.2 

2.0 


-J 



the result of the generalized sign algorithm la 


11.94930 

- 2.24532 

15.31728 

21.65328 

- 2.24532 

- 3.84265 

0.49866 

- 4.59065 

- 7.18665 

0.49866 

- 4.07999 

0.55999 

- 4.91998 

- 7.59998 

0.55999 

- 4.03465 

1.04266 

- 5.59865 

- 7.02664 

1.04266 

4.15732 

- 0.50132 

4.90931 

7.81331 

- 0.50133 


the set of elgenprojectors for the sign spectr'^l decomposition are 



2.08533 

1.04266 

0.52132 

5.21332 

1.04266 


- 0.93866 

- 0.46933 

- 0.23466 

- 2.34666 

- 0.46933 

at 

- 0.95999 

- 0.48000 

- 0.23999 

- 2.39999 

- 0.48000 


- 0.36266 

- 0.18133 

- 0.09066 

- 0.90666 

- 0.18133 


1.06133 

0.53066 

0.26533 

2.65333 

0.53066 


- 9.86397 

3.28798 

- 14.795956 

- 16.43995 

3.28799 


2.90399 

- 0.96799 

4.35598 

4.83998 

- 0.96799 

- 

3.11999 

- 1.03999 

4.67998 

5.19998 

* - 1.03999 


3.67199 

- 1.22399 

5.50798 

6.11998 

- 1.22399 


- 3.09599 

1.03199 

- 4.64398 

- 5.15998 

1.03194 


8.51195 

- 6.46416 

14.40820 

10.95989 

- 5.66391 


- 2.03199 

1.90406 

- 4.08809 

- 2.55997 

1.10396 

= 

- 2.15998 

1.52006 

- 3.44009 

- 2.79996 

1.51996 


- 3.17598 

2.47203 

- 5.48403 

- 4.07996 

2.07198 


1.96798 

- 2.09607 

4.41211 

2.43996 

- 0.89595 
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0.26668 

2.13350 

- 0.13358 

0.26674 

1.33326 

0.06666 

0.53326 

- 0.03323 

0.06664 

0.33336 

> 0.00000 

- 0.00006 

0.00010 

- 0.00002 

0.00003 

- 0.13334 

- 1.06670 

0.06671 

- 0.13335 

- 0.66665 

0.06667 

0.53341 

- 0.03345 

0.06669 

0.33329 


It can be shown this set of elgenprojectors satisfy the properties of the 
I elgenprojectors given earlier. 

I 

I 

I 

I 
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I 
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I 
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3. A Computing the pth Root of a Positive Semidefinite Matrix 


Several recent papers have appeared on the computation of roots of Batrlcaa 
that are generalizations of Newton's method. Hoskins and Walton, [5], described 
a method of computing the pth root of a positive definite matrix with p an in- 
teger. Denman [17] described an extension to the Hoskins and Walton algorithm 
with the only restriction that the matrix not have eigenvalues along the axis. 
An algorithm will be given in this paper that removes the restriction on the 
presence of zero eigenvalues. The procedure can be extended for complete rsmoval 
of the restriction on eigenvalues along the J(<^-axis. 

Let A be a general positive semidefinite nxn matrix with the possibility 
of a zero eigenvalue. Define A as the p-th root of A such that 

(3.4.1) A - A** 

The Newton algorithm as given by Hoskins and Walton is 

/\ 

which will converge to A in the limit as e-H) with stopping criteria 
(3.4.3) ||X,^,-XJ|<C 

where c is a preselected error bound. It can be seen that the above iterative 
algorithm fails when there is a zero eigenvalue. 

The .generalized sign matrix provides the mathematical tool to eliminate 
the problem due to the zero eigenvalue. If A is positive definite, then the 

34 
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sign of A will be I. The generalized sign matrix will give 


(3.4.4) 


S(A) - P 


from which it follows that the null eigenprojector is given by 


(3.4.5) 


0 

P - I-P - I-S(A) 


The null eigenprojector will have the same eigenvectors as A but will have 


eigenvalues -t-1 for all zero eigenvalues of A and zero eigenvalues for all 


R^(A^)>0 of A. Since the eigenvectors of P are the same as those of A, 


the sum of A and will have the same eigenvectors as A but A+P^ will be a 


positive definite matrix with -t-1 eigenvalues in place of the zero eigenvalues 
cf A. If A^ is defined as 


(3.4.6) 


S - A+P' 


then the p-th root of A^^ will be correct aside from the inclusion of the p-th 


root of 1 rather than the p-th root of zero. The root due to P is removed 


1/p. + 

from the p-th root pf A^^ by multiplying A^ by P or 


(3.4.7) A = * (A+p‘^J‘^V 


An example is given wliere A is 4x4 and its eigenvalues are ■ 0, X2 “ 0* 






Xj = 3.258 and X^ 


10.741 
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A 


2.0 - 1.0 1.0 

-1.0 4.0 3.0 

1.0 3.0 4.0 

-1.0 -3.0 -4.0 

.a 


- 1.0 

-3.0 

-4.0 

4.0 


Using p ■ 0.5, Che generalized sign matrix Is obtained with 7 iterations for 


the right and left shifts when e " 

lO"^ 

with 



” 0.6 

-0.4 

0.2 

-0.2 


-0.4 

0.6 

0.2 

-0.2 

8 (A) ■ P ■ 

0.2 

0.2 

0.4 

-0.4 


-0.2 

-0.2 

-0.4 

0.4 


I 

I 

I 

I 

i 

1 

I 

I 

I 

I 

I 


The square root, p-2, was obtained from the Newton algorithm in 5 iterations 
with 





1.091906 

-0.662353 

0.429553 

-0.429553 


-0.662353 

1.485410 

0.823057 

-0.823057 


0.429553 

0.823057 

1.252610 

-1.252610 


-0.429553" 

-0.623057 

-1.252610 

1.252610 


The 5-th root was obtained for the same matrix and 


converged in 7 iterations 




= A ■ 


0.761942 

-0.492761 

0.269181 

-0.269181 


-0.492761 

+0.853145 

0.360384 

-0.360384 


0.269181 

0.360339 

0.629565 

-0.629565 


-0.269181 ” 
-0.360384 
-0.629565 
0.629565 


The error in reconstructing A from A 


1/2 


and A 


1/5 


was no greater than lE-05 
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for all eleiiienCh, Chat la 




I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

a 

I 

I 




I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 
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3.5 ComputaClon of th« Generalized Inverae Matrix of A 

The generalized inverae matrix can alao be obtained from the eigen- 
projectors as given earlier. The null eigenprojector will be used in the 
same, manner as it was used for finding the p-th root. 

Let A be an n><n matrix with zero and nonzero eigenvalues and can be 
represented as 

(3.5.1) A - M{dlagIX.,0,0 A .X 

J. r TTX n — 1 n 


The generalized inverse matrix is obtained when the nonzero eigenvalues are 
replaced by the reciprocal with the zero eigenvalues unchanged. If A^ is 
the generalized Inverse matrix with 


(3.5.2) A^ - M{diag(X"^,0,0,...,0,X”^,A"^^ 


it follows that A satisfies the four properties of Penrose [26] which are 


3.5P.1 


3.5P.2 

3.5P.3 

3.5P.4 


AA^A 


a'*’aa^ 


A 


t 


(A^A)^ - A^A 

+ 1 t 
(AA^)^ = AA^ 


Assume that the generalized sign matrix of A has been found. It then 
follows that 

(3.5.3) A - A+P° 
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Thl0 matrix 


where has been computed as in (3.3.9) for the most general case, 
will have rank n since 


(3.5.4) n ■ rank(A) + rank(P^) 


A 

and is invertible. The inverse of A is given by 


(3.5.5) 


A"^ - (A+P°)“^ - A''‘■^P° 


1 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


Since the eigenvectors of P^ are the same as those of A. The generalized 
inverse of A is then given by 


(3.5.6) A^ - 


'S— 1 

or A - A (I-Pq). 

+ 

The matrix given in Section 4 was used to compute A by this method with 


2.41U7 

19.45559 

-30.77219 

8.67789 

-11.74442 

-1.18890 

- 4.34446 

6.32779 

-3.12226 

2.85556 

-1.00002 

- 5.00001 

7.49999 

-3.00004 

2.99999 

-0.78890 

- 7.14446 

11.42778 

-2.92226 

4.95554 

1.47780 

3.98890 

- 5.50555 

3.54449 

-2.31111 


U i;an bo shown that the eigenvalues of this matrix arc given by • 0, 
= -0.5, = 0.33333,- = 




-J and Xj ■ + j • 


3.6 Decoupling r>£ Linear Time- Invariant System Equations 


The last topic to be discussed in this paper is the use of the spectral 
decomposition concept for solving system equations of the form 

(3.6.1) x(t) - Ax(t) x(0) ■ X eR^ 

where A is nxn and x(t) is nxl. The solution x(t) can be expressed in the 
form 


(3.6.2) x(t) =» <t'(t,0) x(0) “ exp(At)x(0) 

where ())(t,0) is the state transition matrix with 
The general matrix A can be decomposed into 

(3.6.3) A = A^+A +A^+A^ 

where A^, A , A^ and A^ are as defined in Section 3.2. It follows that (3.6.2) 
can be modified with 

(3.6.4) x(t) = exp[(A^+A +A^+A^) t]x(0) 

= exp(A^t)exp(A t)exp(A^t)x(0) 

since A^ = 0. The exponential matrices are commutative as each matrix has the 
same eigenvectors. 

The initial condition vector can also be decomposed with 
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P 


0.07692 

-0.23076 

0.69230 

-2.07692 


-0.07692 

0.23076 

-0.69230 

3.07692 


0.01923 

-0.05769 

0.17307 

-0.51923 


-0.01923 

0.05769 

-0.17307 

0.51923 


P 


1 


0. 32307 
-0.36923 
-1.29230 
-1.29230 


-0.12307 

0.56923 

0.49230 

-2.27692 


-0.16923 

-0.09230 

0.67692 

0,36923 


-0.03076 

-0.10769 

0.12307 

0.43076 


P 


0 


0 


the resultant eigenprojectors on the initial condition vector are 


"0.2*" 


~ 0.307692 


-0.892308" 


0.2 

_ 

0.923077 


2.876923 


0.2 

X (0) = 

-2.769231 

It 

z> 

X 

3.569231 

)( (0) " 

0.2 


8.307693 


-11.507692 



The initial condition vectors given above were used to compute the vectors 
x^(t), X (t) and x^(t). The different modes were separated by the algorithm 
according to tlie spectral domains. Table 3.1 gives the third element of 
Iho vectors ;is well, as the third clenHint of x(t) as computed directly from 


(3.6.1) 
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r 


1 

t (sec) 

X3(t) 

X.3(t) 

1 

X 

xjct) 

1 

1.000000 

0.200000 

-2.769231 

3.569231 

1 

-2.176561 

0.298365 

-0.834076 

-1.640849 

■ 0.8 

-5.661723 

0.445108 

-0.251218 

-5.855612 

1 

-5.930082 

0.664023 

-0.075665 

-6.518440 

1 

1.6 

-2.259451 

0.990606 

-0.022789 

-3.227268 

1 2.0 

3.492468 

1.477811 

-0.006863 

2.021520 

2.4 

1 

8.246650 

2.204635 

-0.002066 

6.044082 

1 2.8 

9.688690 

3.288929 

-0.000620 

6.400383 

1 

7.780616 

4.906506 

-0.000183 

2.874297 

■ 

3.6 

4.924292 

7.319647 

-0.000049 

-2.395300 

1 4.0 

1 

1 

4.707674 

10.919630 

-0.000006 

-6.211941 


■ All computations were in single precision. The value of x^(t) should 

+ - I 

I be equal to tlie sum of x^(t) , x^(t) and x^(t). The difference between x(t) 

and the sum of the individual terms was 9x10 ^ at t =• 4 seconds. 

I 

I 

I 

I 

I 

I 
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4. Laplace Transform of Matrix Functions 

Assume that a system Is defined by a set of differential equations of 
the form 


(4.0.1) x(t)+Aj^x(t)+A2x(t) ■ Gu(t) 


where A^, A^ and G are constant matrices. The matrices A^ and A^ are nxn 
and G Is nxm. Since this type of equation Is encountered In structures, 
the closed form solution In the time domain Is of Interest. The Laplace 
and Inverse-Laplace transforms are therefore of Interest. 


4.1 Laplace Transform 

Consider a system defined as In (4.0.1) with Initial conditions x(0) 
and x(0) from which It follows that 


(4.1.1) 




r 

M 

O 

1 


x(t)’ 


1 

o 


■ 




+ 


x(t) 


.■■^2 ■■^1 


x(t)^ 


LgJ 


u(t) *Ax(t)+Gu(t) 


with 


sX(s)-x(0) 

s^X(s)-sx(0)-x(0) 


."^2 "^1. 


X(s) 

sX(s)-x(0) 


0 

G 


u(s) 


The second equation of (4.1.2) gives 

(4.1.3) s^X(s)-3x(0)-x(0) - -A2 X(s)-Aj^sX(s)+Aj^x(0)+GU(8) 
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which can be rearranged with 


(4.1.4) [8^I+Aj^8+A2]X(a) = t sl+Aj^]x(0)+x(0)+GU(s) 

2 -1 

Premultiplying by [s I+Aj^8+A2] givee 

(4.1.5) X(e) “ [s^I+Aj^8+A2] [sI+Aj^]x(0)+x(0)+Gu(8) } 


Suppose that the Inverse of [sl-a] Is now considered where 



• 

-1 

“ 2 . 

a 

-1 

“ 

■ 

si 

-I 


8 I4-Aj^s+A2 

0 


sI+Aj^ 

I 

A- 

sI+A, 


0 

8^I+A,8+A« 


-Ao 

si _ 

L 2 

IJ 


• 

1 2 J 


^ 2 



which will be written In the compact form 


(4.1.7) [sl-A]'^ . [s^I+Aj^ 8 +A 2 ] 

for convenience. It should be noted that the upper row of [sl-a] ^ contains 
the two terms of Interest In (4.1.5). 

If (4.0.1) has a leading term A^, then (4.1.5) becomes 


sI4Aj^ I 

-A2 SI 


(4.1.8) X(s) - IAq 8 ^+Aj^ 8+A2 ]"^{[SI+Aj^]AqX( 0 )+AqX( 0 )+Gu(s)} 


which requires that [si-a] be modified. A general form for A^ In [sl-a] Is 
not known thus (4.1.7) will be used with A^ Included as shown In (4.1.8). 


4.2 Inverse Laplace Transform 

The Inverse Laplace of the desired matrix equation can be obtained 
from (4.1.7). Consider the eigenvalue-eigenvector problem for A where 


(4.2.1) 


XI- A 


XI -I 

Xl+Aj^ 


where X and s are Interchangeable In the mathematical operations. It Is 
known that the eigenvector matrix of (4.0.1) has the general form 


(4.2.2) 






d> A 

^2 2 . 


with 4>2 " and A2 * A^* for a system with all eigenvalues In pairs, real 
or complex conjugate. Since the structure problem Is of that type, then 9 
can be assumed to be 


(4.2.3) 


4 - 




* 

4>1 

^11 


It then follows that 


(4.2.4) 


■ *1 


> 

(-• 

0 


r * 1 

.♦1*1 ♦x. 


1 

0 

> 

M » 


. .*,* 


Using the property that a function of a matrix Is given by 
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(A. 2. 5) f(A) - 

I 

I 

then 

I 

I (4.2.6) exp(At) - ~^[sl-a] - W(s)4>"^ 


I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


from which It follows that 


(4.2.7) exp (At) ■ $ 


exp(Aj^t) 


exp (At) 


.-1 


It Is not difficult to show from (4.1.8), (4.2.4) and (4.2.7) that 


(^•2-8) X ^Is^I+Aj^ 8+A2] ^ - exp(Bj^t)sin(B2t)B“^ 

(4.2.9) ■35^^[s^I+Aj^s+A 2]"^ - exp(Bj^t)[cos(B2t)-sln(B2t)(B2^^)] 

2 

• • ® spectral factor of [s I+Aj^s+A 2 ] or 

(4.2.10) [8l+Bj^+jB2][SI+Bj^-JB2] - [s^I+Aj^s+A2] 

with upper half plane and Bj^-jB 2 In the lower half plane. 

Substituting back into (4.1.5) with A^ - I, x(t) is given by 

(4.2.11) x(t) - {exp(B^t)[cos(B2t)-sln(B2t)(B”^Bj^)]}x(0) 

+ (exp(Bj^t)sln(B 2 t) (B“^}x(0) 

+ j exp[Bj^(t-T)]sln(B2(t-T)](B2^)ru(T)dT 





"Page missing from avaiiabie version" 


4.3 Facorlzatlon of A(s) 


Let A(s) be defined as 

A(s) “ 1 s^+Aj^8+A2 

for which the factors [sI+B,+JB 2 ] and [sI+Bj^-jB^] are to be obtained. It 
can be shown that there exists a matrix R such that 

(4.3.1) R^+A^R+A2 - 0 
where R Is given by either 

(4.3.2) R^ - - -B^-jB2 


or 


(4.3.3) R 2 - - -Bj^+JB2 


The matrix R will be complex and can be computed from the eigenvalues and 
eigenvectors of A where 


(4.3.4) 



I 



It can also be shown that A(s) can be written as 
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(4.3.5) 


A. 2 A 

A(s) - (is+ -^y ^ 


which factors as 


(4.3.6) A(s) - {l8+ ^ +JIA2 -aJ/4]^^^}{Is+ ^ -j [A2 -aJ/4]^^^} 


therefore 


(4.3.7) 




(4.3.8) 


B 2 - IA2-aJ/4]^^^ 


provided that A2~A^/4 Is a positive definite matrix. Provided that this 

re itrictlon Is satisfied, It Is not necessary to compute the algebraic 

2 

solution to (4.3.1). The square root of k^k^/l* can be computed by the 
procedure given by Denman 117] or Hoskins and Walton [20]. 

The equation given In (4.3.1) Is a algebraic matrix Rlccati equation 
and plays an liiq>ortant role In system analysis and particularly In spectral 
factorization. To show that R Is as given (4.3.2) consider (4.3.1) and 
(4.3.4). It Is kno%m that the eigenvector of a Is as given In (4.2.2) thus 


(4.3.9) 


T r . *1 

■I 

[+aJ 


thus If X > then 


(4.3.10) A2<|>j^+Aj^<(»j^Aj^4Aj^<j)j^A^ - 0 
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Now if 

(4.3.11) 
therefore R 

(4.3.12) 

thus 

(4.3.13) 

(4.3.14) 
where it is 

(4.3.15) 

If • 

(4.3.16) 


• then 

(|.iAJ<|»‘^+Ai(()iAi4)"Va 2 - 0 
- as given in (3.2). Now from (2.4) 


' 0 I ■ 


■ <^1 K ' 

A “■ 


■ h 

A * •* 

♦1 


“ 

A A 

.* 



A A 

-A2 -Aj^_ 


1^*1 

0 A, 
IJ 


• Vi 



Ai - -(4 .^Aj4»‘^-(|)*aJ‘(|)J"^)(<|)^A^4.’^-<|.JaJ((>J"^)'^ 

necessary to prove that 

Rl - - +J(A2 "aJ/ 4)^^^ - ♦l^l'*’!^ 

•k 

R2 “ -Bj^+jB2 then 

Ai - -[(-B^+jB2)^-(-B^-jB2)^][-B^+jB2-Bj^+jB2]‘^ 

- -(Bj^B 2+B2B^(B2)"^ - -2B, iff B^B2“B2B^ 

A2 - -[(-B^+JB2)-(-B^-jB2)][(-B^+JB2)"^-(-B^-JB2)"^l‘^ 

- -2jB2(-B^-jB2)(-2jB2)‘^(-Bj+jB2) - B^+B^ iff BxV®2®l 


(4.3.17) 


I 
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Substituting back into for and A^. 

2Bi 222 

( 4 . 3 . 18 ) 


where 


( 4 . 3 . 19 ) 



®2 “ 


I 

I 

! 

I 

I 

I 

I 

I 

I 

I 

I 

I 


5. Olagonalizatlon of Matrix Polynomiala 


There haa been niimerous conments made at technical meetings regarding 
the dlagonallzatlon of matrix polynomial of the form 

(5.1) A(s) - M3^-H:8+K - M^^^[l8^+Cs+K]M^'^^. 

It haa been ateted that the aecond order matrix polynomiala of the 
general form cannot be diagonalized Into the form 

(5.2) D(8) - 

where 2 qu ) and u) are diagonal forms. It Is true that there does not exist a 
constant matrix such that 

(5.3) D(s) QA(s)Q"^ 

but It is not true that A(s) cannot be diagonalized. As will be shown In 
the following sections, polynomials do exist such that 

(5.4) D(s) Q(8 )A(s)P(s) 

To show tl'ls. It will be assumed that the mass matrix Is Invertible so that 

1/2 

M exists. This assumption may be too restrictive but will be used. 

5.1 Dlagonallzatlon of Block Companion Matrix 
Since it has been assumed that M Is Invertible, a block companion form 


matrix a can be defined with 



■^3 


(5.1.1) A - _ 

l-K -C 

where K - •nd C - If K and C are ayimatrlc then K 

and C will be symmetric when Is symmetric. 

The matrix A can be de/’lned as 

(5.1.2) A - 


where ^ Is the n^n matrix of eigenvectors of a and A Is the eigenvalue 
matrix. Let A have distinct eigenvalues with real and complex eigenvalues 
In conjugate pairs. The eigenvalues and eigenvectors of A can always be 
ordered such that 


(5.1.3) 




If “k k 

dlagt^l^l^2^2 ' • * ^m^m^mfl^mf2 ' ‘ * ^2n^ 


. r * 1-1 


'2n^ 


I 


and thus 


(5.1.4) 


A » *”^A4 - dlag[Aj^A*A2^2” *^2n^ 


Consider any pair of eigenvalues Including the real eigenvalues. The com- 
parison form of the 2x2 subblocks can then be written as 


(5.1.5) 


• ‘ 1 • \l 

7*^1 "^^i“iJ L^i ^1 






-1 





54 


with a new matrix defined as given by 


FA 


(5.1.6) Ag^ . 


B1 


B2 


"Bill 


hco. 




0 

0 

0 

2 

-U'o 


0 

0 

1 

■^^2“2 


0 

0 

0 

0 


-U) 


n 


0 

0 

0 

0 


-2C u) 

n n 


It follows that a row-column Interchange can now be made to t ansform 
Ag^ into the form 


(5.1.7) 


DC 


— U) —2 Qii 




where E Is an elementary matrix. For exaiq>le, consider the 4x4 matrix A 


with 


BC 


(5.1.8) 


B1 


0 

0 


■2^l“l 

0 

0 


0 

0 

0 

2 

-Wo 


0 

0 

1 


- 2 ^ 2“2 


then 




10 0 0 
0 0 10 
0 10 0 
0 0 0 1 


-0.J 

0 

0 


-2C 

0 

0 


0 

0 

0 

-J. 


0 

0 

1 




10 0 0 
0 0 10 
0 10 0 
0 0 0 1 


(5.1.9) 
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The transformation process requires three transformations with 


( 5 . 1 . 10 ) 


where E Is the elementary matrix, G Is a block diagonal matrix with 


0 


1 




or a suitable pair of real eigenvalues and ^ Is the usual eigenvector matrix. 
Yh'.'! 0^ matrix must be Invertible which Is a valid assumption for conjugate 
pairs of eigenvalues and distinct real eigenvalues. 


5.2 Matrix Polynomial Form 


The analysis given in Section 5.1 was concerned with the "state vari- 
able" or block companion form and does not give the required transformation 
in matrix polynomial form. The objective of the dlagonallzation form is 
to obtain Q(s) and P(s) as given in (5.4). The desired Q(s) and P(s) can 
be found from the work in Section 5.1. It is not difficult to show that 


(5.2.1) 


[A(s)] 


03 

M 

1 

M 

1 

-1 

1 

iH 

m 

\< 


• • 

si-fC 1 

L K SI4C 


0 [A(s)]"^ 

* sJ 


-K si 

m q 


which was shown in a previous report . The dlagonallzation process given In 
the previous section holds for a(s) with 


(5.2.2) 


E04>’^[A(s)]*0"V 


si 
2 

oj s 1 + 2 ^ 0 ) 


■1 

2C03J 


with 


(5.2.3) 


81 

2 

0) 


-I 
8I+2CU) 


T -1 


[D(s)]’^ 


[D(s)] 


-1 


sl+2?u) 

-0,2 


I 

si 


E0»‘^[A(s)]'^^"V 


Suppose that the product E0$ ^ Is defined as Q ^ thus 


(5.2.4) 


[D(s)f^ 


• « 

sl+2^0, I 

«“1 

iA(s)]"^ 0 


1 — 
M 

M 

0 [D(s)]"^^ 


-0,2 8 1 

m Ml 

- Q 

. 0 [A(s)]’^ 

• mi 


1 

1 

cc 

M 
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(5.2.5) D(a> - [(8l+C)Q^2'^22^ "[A(*) ] [sQ^2'^^l^ 


as well as 


-Iv. 


(5.2.6) D(s) - [8Q22-*^Qi2^ A(s) [sQ22+<J2iJ 


for the quadratic matrix polynomial with 


(5.2.7) 


Q - 


^11 ^]2 


L ^21 ^22 


where Is an n^n partitioned block of Q. It follows from further analy- 
sis that 


(5.2.8) 


0 )^ - "^ 12^21 “ ^‘^^ 12 ‘^ 22 ^ ^^11 


and 


- 1 , 


(5.2.9) 2CW - Q22lQ2l‘*^^22"‘“^12^ 


" ^ 12 -^ 22 "^ 11 ^ 


The usual transformation to normal coordinates 


(5.2.10) x(t) - *q(t) 


Is replaced by a coordinate transformation of the form 
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(5.2.11) Q22q(t)+<)2j^q(t) - x(t) 
thus 

(5.2.12) [l8^+CsfK][Q228+Q2i]q(s) - U(s) 

If the force vector u(t) Is transformed with 

(5.2.13) Q22ht)-KQ^2^<‘'> “ 
then the resulting equation Is 

(5.2.14) [l8^+CsfK][Q228-K)2^]q(8) “ lQ228-KQj^2^^^®^ 
therefore 

(5.2.15) [Q228-i^j^2^"^f^®^'^®^'^^f‘^22®‘^21^‘^^®^ " ^^®^ 

It follows from (5.2.6) that 

(5.2.16) D(s)q(s) - f(s) ■ [ls^+2Cws-Ho^]q(s) 

A different form of (5.2.15) can be obtained from (5.2.5) although 
Is the same. 

The mass matrix must now be taken back Into the equation since 


5.2.16) 

was 


factored out. 


5 . 3 Examples 


As an example of Che above procedure, let 


A - 


0 

0 

1 

0 

0 

0 

0 

1 

-35 

-40 

12 

8 

-40 

-35 

8 

12 . 


with A ■ [15 5 5 -1] and 


-1 

-5 

103 

-1 

-1 

-4 

-23 

1 

-15 

-25 

515 

1 

-15 

-20 

-115 

-1 


as conq>uted with by EISPAK with column vector scaling. The transformed 
matrix was conq>uted as 


with 


and 



-1.9E-6 

-2.30E-6 

-75.0 

2.38E-7 


-2.86E-6 1 -4.76E-6 

2.62E-6 2.98E-8 1 

-1.22E-4 20. -7.63E-5 

5 -7.45E-9 4 



^ 12^21 


-75. C -3.81E-6 
-1.19E-7 5 
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201 - QI^022-‘*11> 


-2.0 3.81E-6 

-2.98E-8 -A 

» 


All computations were In single precision except for the eigenvalue. The Q 
matrix required for the block dlagonallzatlon of A was computed and Is given 
by 


-7. 

16.3333 

O.A 

17.333 

-5.5 

- 3 

0.3 

- A.O 

-30 

86.6667 

1 

85.6667 

-22.5 

-20 

0.5 

-19 


with 

Q - 4>0"^E . 


The example Indicates that the block dlagonallzatlon of A Is valid. 
The matrices 9,Q, and E will alw&ys exist thus Q will transform A Into a 
diagonal block companion matrix. 



6. Quadrature Method and Laplace Transforms 

The quadrature formula for numerical Integration Is a useful mathe- 
matical tool for system analysis. The primary use In the work described 
In this report Is that related to system Identification although the 
procedure has other related uses. Readers that are not famllar with the 
work should consult Bellman, Kalaba and Lockett for details as well as 
applications, [27]. Numerous papers are also available In the literature 
on the use of other orthogonal polynomial rather than the Legendre poly- 
nomials as used In [28]. 

Consider a function f(t) that Is Laplace transformable, that Is 


( 6 . 1 ) 


r 


I f (t) !e ^^dt<° 


ftr a real positive o. The Laplace transform of f(t) Is given by 


( 6 . 2 ) 


r 

■*0 


F(s) ■ i f(t)e ®*"dt 


where the Inverse Laplace transform Is given by 


(6.3) 


f(t) 


2tt1 


fO+±Ui 


a-lo) 


F(s)e®*^ds 


The Laplace transform of f(t) will be denoted by F(s) with the notation 


(6.4) gt[f(t)] - F(s) 


and the Inverse transform by 
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(6.5) X'^P(b)] - f(t) 

The Laplace transform and the Inverse Laplace transforms also hold for 
vector and matrix functions. If f(t) Is a nxl vector then F(s) will be a 
nxl vector where 

(6.6) F^(s) - f^(t)e"®‘'dt 

with a similar expression for the Inverse transform. Let A(t) represent a 
matrix, then the elements of the matrix have Laplace transforms A^^(s) de- 
fined by 

(6.7) A^j(s) - I” A^j(t)e'®‘=dt 

The operations are always performed on each element of the matrix rather 
than as a matrix operation. 

6.1 Properties of Laplace Transforms 

Some liq>ortant properties that will be useful In the work that follows 
will be given 

a) Real differentiation. If ;{[f(t)] ■ F(s) then 
(6.1.1) X iffl - SF(s)-f (0-) 


and 
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I 

I 

I 

E 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


(6.1.2) - ■"F(8)-8"‘^f(0-)-8"""h0-)... 

dt" 

where f(O-) indicates the first derivative with respect to time and 
indicates the (n-l)th derivative. The argument 0- represents the initial 
condition on f(t). 


b) Real integration. Let '^[f(t)] - F(s) then 


(6.1.3) 'X.l 


ft 


f(T)dl] 


. lisi 


c) Differentiation by s. If s is the Laplace variable, then 


(6.1.4) (tf(t)] - - 


6.2 Laplace Transforms and Quadrature Integration 

Let f(t) represent a Laplace transformable variable, scalar or vector, 
then [f(t)] ■ F(s) with 

(6.2.1) F(s) - f f(t)e"®*"dt - I w.r®“^f(t.) s - 1,2,3 « 

Jo 1-1 ^ ^ ^ 

where the summation form is the quadrature integration of f(t). The weights 
are denoted by w^ and the roots by r^ with t^ - -log(r^). Details of the 
derivation of (6.21) are given in [27]. The weights and roots are given in 
the appendix to this report for N - 3, 4,..., 15. Extensive tables can be 
found elsewhere [29]. 
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The properties given in Section 6.1 can be used to derive 
other ' seful expressions. From property a) 


( 6 . 2 . 2 ) 


? .-1 

»F(.) . f(0 ) - I w r’ 
1-1 


with t(0 ) the Initial condition on the function f(t). From b) 


(6.2.3) 


lisl 

8 


r 8-1 , 

I w^r [ 
1-1 ^ ^ 


ft. 


f(T)dx] 


and finally 


(6.2.4) 


ms) 

ds 


N 


1-1 ^ ^ ^ ^ 


As examples of each of these, let 


(6.2.5) 


f(t) - 1-e 


-t 


which has the Laplace transform of 


(6.2.6) F(s) - 


8 (s+l) 


with 


(6.2.7a) 


ft 


(6.2.7b) 


(l-e”^)dT - t+e“^-l 


(6.2.7c) 
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tf(t) - 

The table below gives computed values from (6.2.1), (6.2.2), (6.2.3) and 
(6.2.4) with exact values for quadrature orders N-7 and N-11. The functions 
F(s) and sF(s) are correct to six digits whereas dF(s)/ds and F(s)/s have 
significant errors. The errors In the latter two functions decrease with 
the quadrature orders. It Is therefore necessary to select the quadrature 
order higher If f(t) Is not a smooth fimctlon over the range of time used 
In the procedure and partlcularlly If f(t) has variations between the 
sampling points. Since the Interval time between the sample times t^ In- 
creases as 1 get larger, less Information Is available for lerge t^. 


s 


F(s) 

sF(s) 

dF(s)ds 

F(s)/s 

1 

EXACT 

0.5 

0.5 

-0.75 

0.5 

(N-7) 

COMPUTED 

0.5 

0.5 

-0.738658 

0.488739 

(N-11) 

COMPUTED 

0.5 

0.5 

-0.745205 

0.495219 

2 

EXACT 

0.166666 

0.333333 

-0.138888 

0.083333 

(N-7) 

COMPUTED 

0.166666 

0.333333 

-0.138972 

0.0834144 

3 

EXACT 

0.083333 

0.250000 

-0.048611 

0.0277777 

(N-7) 

COMPUTED 

0.083333 

0.250000 

-0.048609 

0.0277758 

(N-11) 

COMPUTED 

0.083333 

0.250000 

-0.048611 

0.0277776 

4 

EXACT 

0.050000 

0.200000 

-0.022500 

0.012500 

(N-7) 

COMPUTED 

0.050000 

0.200000 

-0.022500 

0.012500 

(N-11) 

COMPUTED 

0.050000 

0.200000 

-0.022500 

0.012500 

5 

EXACT 

0.033333 

0.166666 

-0.012222 

0.006666 

(N-7) 

COMPUTED 

0.033333 

0.166666 

-0.012222 

0.006666 

(N-11) 

COMPUTED 

0.033333 

0.166666 

-0.012222 

0.006666 


Table 6.1 COMPUTED FUNCTIONS 


H The quadrature formula given In (6.2.1) falls for noninteger values of 

_ s. It has been shown In [27] that values of F(s) other than the Integer 

I 

I 
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valuaa can be coaputed by aodifylng the formula. Let a be an arbitrary 
scalar , then 

N . 

(6.2.8) F(a/a) - a ] w r“'^£(at.) 

1-1 ^ ^ ^ 

Equatic«<e (6.2.1) and (6.2.8) can be used to compute F(s) such that rapidly 
changing values of f(t) are included in the data set for F(s). This pro- 
cedure will assure that variations of f(t) for large t^ has been included 
in the s-domain information. 
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6.3 Inverse Lsplace IransfomiB and Quadrature Formulas 

The time function f(t) at times t^ can be computed from F(a) by noting 
that a set of linear equations can be defined from F(s) with 

(6.3.1) F(s) - y w r®”^f(t.) 

1-1 ^ ^ 1 

Let S - 1,2,3,. ..,N, then 




’f(i)' 


■ 

w 

1 

''2 

"3 

• • • 

m 

w 

N 


' f(t^) ” 

(6.3.2) 


F(2) 

■I 

w r 
11 

• 

” 2*^2 

• 

” 3*^3 

• • • 

w r 
N N 


f(t2) 





N-l 

w r 
*11 

• 

N- 

” 2*^2 

■1 N-l 

”3'^3 

• • • 
• • • 

N-l 

"n*^n . 



Equation 

(6.3.2) can then 

be considered 

as having the 

form Ax - b 

which can 

be solved by Inverting A with x - 

A~\. 

The A matrix 

is 

ill-condlti ~>ned and 


will require care in selecting the algorithm to conq>ute f(t). 

The computed values from (6.3.2) and the exact values of the f(t) given 
in (6.2.5) are tabulated in the table below. Single precision calculations 
were used with a 7th order quadrature formula. The singular value decom- 
position was used to compute the f(t.) values. 


1 

h 

(sec) 

f(t; 

Exact 

L> 

Computed 

Error 

1 

0.025772 

0.0254458 

0.0254146 

3.12E-5 

2 

0.138382 

0.129234 

0.129285 

-5.1E-5 

3 

0.352509 

0.297077 

0.297007 

7.E-5 

4 

0.693147 

0.5 

0.50009 

-9.E-5 

5 

1.213762 

0.702923 

0.702812 

l.llE-4 

6 

2.046127 

0.870766 

0.870889 

-1.23E-4 

7 

3.671195 

0.974554 

0.97445 

1.04E-4 


Table 6.2 Compute Value of f(t^) from F(s) 
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Additional values of f(t^) could be obtained by using a higher order quad- 
rature <'<}rmula or by utilizing (6.2.8). The accuracy of the computed values 
of f(t^) will depend upon the order cf the quadrature formula as well a^ 
the variation of f(t). There Is no assurance that a low order quadratu^ 
formula will give a unique vector f(t) when the s-domaln Information Is In- 
adequate. For example, a resonance peak may be missed entirely or else the 
8-domaln Information may not cover a region In which F(s) varies rapidly as 
at cutoff. The same Is true for the computation of F(s) from f(t). 

The above procedure can be utilized to conqiute the response of the 
structure to a forcing fun tlon. Let X(s) be defined by 

(6.3.3) x(s) “ [Ms^+CS+K] ■^Gu(s) 

2 -1 

where [Ms -t-CS-t-K] can be expanded as 
2 —1 

(6.3.4) [Ms +CS+K] - L - 7 ^ R, nxn 

1=1 s+s^ 1 

where Is the matrix residue of the root s - -s^. The residues have .rank 
1 and arc given by the outer product 

9 

(6.3.5) R^ - y^z^ 


/\ 2 '' 

where y^ la the right latent vector of Ms +CS+K and z^ Is the left latent 

2 /N /S /V ^ 

vector. If MB +CS+K la syametrlc then - y^. The vectors y^^ and z^^ust 
be properly normalized. Arsisnlng that the roots s^ appear In conjugate 

* •k 

parrs, the residue for s^ must be given which means that all of the 



2 

residues csn be stored In n locations, n vectors of dimension nxl. An 
additional storage of n locations must be available for the n roots. 

The values of X(s/a)/a can therefore be confuted by using 

2n R 

(6.3.6) X(s/a)/a - ^ • Gu(s) 

It follows that x(t^) can be computed by using (6.3.2) with 

N 

(6.3.7) X(-alnr ) - I P X(s/a)/a - X(t ) 

^ 1-1 ^ ^ 

where Is the first row of the Inverse of the matrix In (6.3.2). Appr&l- 
mately n^4n storage locations will be required by using the above algorism 
which exceeds the storage requirement for a 4th order Runge-Kutta when the 
sparse matrices M, C, _..d K are handled properly in the Integration algorithm. 
The sparse nature of M, C and K cannot be used In the quadrature algorithm 



as described. 



7. Identification of Linear Systems by Quadrature Method 


The identification or parameter estimation problem in linear system 
theory is that of determining the system parameters from a sequence of 
input-output data. In some cases, only statistical prop' *'ies of the input 
are known with measured output data available. It will be ■ lumed that the 
inp 'i(t), and the output, y(t), are known. Since the identification prob- 
lem 3f interest is to estimate the parameters of a vibrating system or a 
structure, u(t) will be the excitation of the system with y(t) measured 
displacements of points or nodes of the structure. The variables u(t) and 
y(t) may be vectors and x(t) will be the vector of all nodes of the structure. 

It will be assumed that the excitation u(t) is that signal applied to 
a set of transducers that apply mechanical forces to the structure whereas 
y(t) will be the output of sensors that measure the displacement of the 
nodes. With those assumptions, the vector x(t) can be defined by 

(7.1) Mx+Cx+Kx - Gu(t) 

where M is the mass matrix, C the damping matrix, K the stiffness matrix 
and G the input matrix. The Zui^tcices, M, C and K will be nxn, G is nxm, 
x(t) is nxl and u(t) is mxl. The oi^tput vector y(t) will be given by 

(7.2) y(t) - Hx(t) 

with H a J-xn matrix and y(t) a i,xi vector. G and H will be considered as 
constant, matrices with known values. 
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The Input and output variables, u(t) and y(t), may be sampled with 

values taken at a fixed sampling rate, £^, or at specified time Intervals 

which are not equal. If f Is fixed, then samples are available every T 

8 

seconds. The nonequal Interval sampling strategy will be to take saiiq>le 8 
at a set of specified times, t^, where 4 T for all 1. The two 

saiiq>llng strategies will be denoted as uniform and nonuniform sampling 
rates . 

If uniform sampling Is chosen, the mathematical theory fcr the linear 
system should be formulated In the z-domaln. It can then be shown that 

(7.3) Y(z) ■ T(z)U(z) “ Hx(z) 

where T(z) Is the transfer function In the z-domaln given by 

(7.4) T(z) - Z{H[M 6 ^-K:s+KJ"^G} - (I-exp(Bj^T) fco 8 ($ 2 T) 

- 3j^32^8in(32T;]z+exp(6j^T) texp(3j^T)-(cos(32T)-3j^32^8ln(3j^T))]}. 
• {z^I-2zexp(3j^T)cos(32T)+exp(23j^T) } . 

The details for the derlvaticn of (7.4) are given In the appendix. The 
major problems associated with this approach Is that the parameters of the 
system are dependent on matrix fimctlons, l.e. exp(3^T), cos ( 321 ^) and 
sln( 32 '^)* Small errors In the elements of the matrix elements may be mag- 
nified In the recovery of 3j^ and 32 * 

Although the nonuniform sampling strategy appears to be less desirable 
than the uniform sampling, there exist unique sampling times that leaves 


I 

I 

I 

I 

I 

I 

I 

I 

I 

II 
I 
I 
I 
I 
I 
I 
1 
I 
I 
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the system In the s-domsln. The system parameters can be estimated by 
application of linear algebra theory. In the Laplace domain, the system 
equations are 

(7.5) y(s) - T(s)U(s) - Hx(s) 
with 

(7.6) T(s) - H[Ms^+Cs+K]~^G 

The quadrature formulas hold If the sampling times are consistent with the 
quadrature Integration procedures. 




7.1 Sampled Values and the s-domaln 


Assume that the Inputs u(t^) and outputs y(t^) have been measured at 
the times t^ specified by the quadrature order N. The system equations can 
then be defined by 


(7.1.1) 


N , 

I y(t ) 

1-1 ^ ^ 1 


(7.1.2) U(s) - I w,r®'^u(t.) 

1-1 ^ ^ ^ 


(7.1.3) Y(s) - I w r®"^Hx(t.) 

1-1 ^ ^ 

where H is assumed to be known. The vectors Y(s) and U(s) can be computed 

from (7.1.1) and (7.1.2) when the sequence of Inputs u(t^) and outputs y(t^) 

are known. The state vector x(t^) can then be computed by use of (7.1.3) 

although not all x(t^) can be determined when H Is ilxm with Y(s) i f[x^(s), 

Xo(s),...,x (s)]. Since y . (s)ouc. (s) , then Y(s) will be Incomplete In the 
/ n 11 

sense that not all values of x^(s) are measured. 

If H Is an Identity matrix, then X(t^) can be determined from (7.1.2) 
with X(s) then equal to Y(s). It follows that 

(7.1.4) [Ms^+Cs+K]X(s) = Gu(s) 
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x’^(l) x’^(l) X^(l) 


u^CDg"^ 

4X^(2) 2X’^(2) x'*^(2) 


T 


U^(2)G^ 



M 



• a • 




• 



C 

■ 


• • • 


,.T 


• 



K 



• • • 




. 

2 T T T 

q X^q) qx\q) x\q) 


T T 

_ u\q)G" _ 


If M, C and K are nxn, then q must be equal to or greater than 3n. Exact 
values of X(s) and X(t^) requires that only 3n rows be used as additional 
Information would be redundant. This Ideal situation will not occur as the 
best measuring device can only resolve the measured values to a certain 
accuracy. Numerical errors will also be Introduced. To overcome these 
sources of errors, a set of overdetermined equations should be used which 
means that q>3n. Algorithms, such as SVD, are then used to compute the 
parameters of M, C and K. A version of the SVD algorithm given In Clarebout, 
[30] has been used with good result, a listing of that program Is given In 
the Appendix. 

Measurements of all of the node displacements Is uneconomical as well 
as unfeasible for a large structure. This means that H will be with 
only H nonzero elements In H under the assisnptlon that % arbitrary node 
displacements are measured. Suppose that the first I nodes are measured, 
thus the first i. elements of x(t^) can be computed. The resulting equations 
from (7.1.3) would then be 



9 m 

Y(l) 


/N 

Wj^H 

W 2 H 

ys 

WjH 



r» A *1 

X(tj^) 

(7.1.6) 

Y(2) 

• 

• 

at 

Vi® 

• 

• 

/V 

w„r,H 

d. 1 

• 

A 

w^r^H 

• • • • 

• a • • 

Vs® 

a 

a 


X(t2) 


• 

.Y(q) . 


• 

N-lr; 
t ''l*’l " 

• 

N-1^ 

W2r^ « 

a a a a 

N-l^i 

w^r^ n • . • 

a 

N -in 
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where 



1 

•H 

X 


•'ll 




^22 ° 

(7.1.7) x(t) - 

X3(tj^) 

• 

/N 

H - 

*^23 


• 


0 

ET 

1 


with the last n-1 columns of H deleted since the elements are zeros. The 

ys 

resulting matrix In (7.1.6) will be q^ln where q>iln. The matrix H Is 
diagonal with sensor gains along the diagonal. As stated earlier, the 
displacement sensors and force transducers are collocated. 

Consider now the mass, damping and stiffness matrices In the form 


(7.1.8) 


M - dlag<m^^ m ^2 ^33 ... 


(7.1.9) C - trldlag<^ 


k 


(7.1.10) K - trldlag 

\ ^ 


12 

C23 . . . . . 

‘’22 

*^33 

21 

*^32 

12 

*^23 

k.^„ 

k__ 

22 

33 

21 

*^32 



nn 



nn 


The first two equations of MS +CS+K would then be 




(7.1.11a) 


(7.1.11b) 


m 228 V2 (8)+C2j^syj^ (s)+C 228 y 2 (8)-C238y^(s;-k2j^yj^ ^“^"•^22^2 

+ k23y3(s) - U2(s) 

where It has been assumed that y^(s) ■ x^(s). 

Equation (7.11a) could then be used to write a set of algebraic equa- 
tions with s as the varying parameter with 


r 


Yid) 


(7.1.12) 


q yj(q) 


y^d) 


qyj^(q) 


y£d) 


4y (2) 2y^(2) 2y^(2) 


y^d) 

yj(2) 


y2d) 

^2 


qy2<q) y^Cq) y2^q) 


■1 

- — 


r n 


“ll 


Ujd) 


‘^11 


u,(2) 


*^12 




Ic 


• 


11 


• 


.*^12 _ 


1 

0 

1 


The values of m^^ will be the same, as will c^^ and for the simple 

structure with evenly spaced nodes, slmllarlly c^^ “ *^J1 ^IJ * ^Jl* 

The above equation will be the only one required for conq>lete determination 
of the structure parameters. 

Suppose now that a beam has equally spaced nodes except for the end 
ones which are located at one-half spacing from the ends of the beam. The 

elements Cj^2* ®21’ *^12 ^21 ®q“®^ to ®22* 

^22* ^22' ^23* ^32* ^23 ^32* first two equations are then used to 

conq>ute the parameters with (7.1.12) being solved first followed by the 
solution of the algebraic equations of (7.1.11b). 

If all of the coefficients In M, C and K are different, a full set of 
measurements will be required. The Important property of the parameter 
identification procedure Is that the structures symmetry can be utilized to 
reduce the dimensions of the matrices thus permitting large space structures 


to be Identified. 



7.2 Numerical Results 
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I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 
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The Identification algorithm described In the previous section was used 
to identify the coefficients of the M, C and K matrices as given In Table 
7.1. The algorithm was first used to identify all 36 coefflcl<>nts by equa- 
tion (7.1.5). The results were fair in the sense that the accuracy was 

4 

approximately one part in 10 or an error of 0.1%. The vectors x(s) and 
u(s) were exact for the numerical test. The resulting matrix equation 
Ax <■ b as given in (7.1.5) was with A a 18x19 matrix, X a 6x6 matrix, and 
U a 6x19 matrix. 

The second test performed used the structure of the matrices M, C and 
K as explained in the text and given by equation (7.1.12). The estimated 
parameters are given in Table 7.2 where the diagonal and tridlagonal ele- 
ments only were computed, all zeros were expressed as zeros. The parameters 
were correct to 12 digits with all computations in double precision on a 
Honeywell 60/60 machine (36 bit word length) . 

The elements of the C matrix were then reduced by two orders of magni- 
tude to test the algorithm for less damping. The results are given in 
Table 7.3. Table 7.4 gives the estimated values when the stiffness matrix 
elements were Increased by one order of magnitude with C the same as given 
in Table 7.4. The estimated values of the parameters in C were in error in 
the eight digit indicating that lightly damped structures will be more 
difficult to identify by the algorithm. 

The test matrices are not representative of any structure. The model 
used was for convenience. Further studies will be made to test the effect 
of noise on the estimated parameters. 

I 
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Table 7.1 The M, C and K Matrices Used in the Identification. 
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Table 7.2 Estimated Values of M, C and K for the Matrices Given in Table 7.1. 


I 
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Table 7. A Estimated Values of M, C and K with Very Light Damping. 


8. Stunmary 


The research effort under NASA Grant NSG-1602 has continued with the 
major objective of developing numerical algorithms for large space structures. 
The research Is directed to finding algorithms for analysis, synthesis and 
Identification of large structures. 

Since the finite element model gives rise to the equation MirfCxfKx^f (t) , 
considerable effort has been expended toward understanding matrix polynomial 
which play an Important role Inthe obtaining the dynamics of systems. In 
addition, the dimensions of the matrix cc<if flclents are lower by one half, 
than the state variable model. Analysis of the structure from the state 
variable model appears to be out of the question when one assumes thst the 
matrices M, C and K, although sparse, may exceed the size of 1000x1000. It 
Is difficult to handle matrices of this dimension and much more difficult 
to manipulate matrices of twice the dimension. 

The decoupling algorithm for analysis has continued but the procedure 
does not appear to have the efficiency required for the desired anal3rsls. 

The research has been redirected with more emphasis on the theory of matr'u 
polynomials . 

This report discusses the decoupling algorithm, the theory of Laplace 
transfonos of matrix polynomials and Identification of the three matrices 
M, C and K. The Identification algorithm Is promising and work will continue 
'?n that aspect of the lar space structure. Data from the NASA LRC beam 
facility will be obtained In the near future and the algorithm will be test- 
ed on this data. Work will continue on the assignment of damping to a 
structure which requires a better understanding of transformations In a 
multl'-dimenalonal space. 
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A.l Derivation of the Discrete Transfer Function 
Let T(s) denote the transfer function In the Laplace domain with 


(A. 1.1) 


T(s) - H[Ms^Cb+K]"^G 

- HM ^^^[ l 3 ^+ A ^ 8 + A 2 l " V^^G 


i^ere A^ and A^ are as defined. The discrete transfer function T(z) is 
given by 


(A. 1.2) 


T(z. - HM^^^(1 -z"^)Z{s"^[Is^+Aj^s+A2]~^}M^'^^G 


The term inside the braces can be expanded Into a partial fraction with 
(A. 1.3) s"^[Ia^+Aj^s+A2]'‘^ = a'^s'^-A'^Es+A^^] [Is^+Aj^s+A2]'^ 

which can be rewritten as 

(A.l. 4) s"^[18^+Aj^8+A2]"^ - a"^{I8"^-(8+Aj^/2)[(Is+Aj^/2)^+A2-aJ/41 

-Aj^/2 [ A2 -aJ/4]’^^^ [A2-aJ/ 4 ] [ (l8+Aj^/2) ^+A2 -aJ/4 

for which the z-trsnjform ie 

A A^ 

-1 7 -1 -1 -1 9 1 I 1/2 

(A. 1.5) Z{s ^[I8 ^+Aj^ 8+A2] } - A2 {z(z-D V[z^I-zexp(- T)co8(A2~ ' T] 

•[z^I-2zexp(- ^ T)cos(A2-Aj/4)^/^T+exp(-A^T)l- ~ [A2-aJ/4]“^^^ 

1/2 2 ^1 1/2 

•[zexp(- -f T)sin(A 2 - -^)^'^T] [z'^I-2zexp(- ^ T)cos(A 2 - -^-) ' T 
•fexp (-Aj^T) ] } 
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Since (z-1) le a scalar, (A. 1.5) can be rearranged* as 

(A. 1.5) Z{s ^[Is^+Aj^ 8+A2] } - A2^{z[z^I-2zexp(- T)cos(A2- 

A A^ 

+ exp(-A^T)]-U^I-zexp(- 2~ T)cos(A2" ^)^^^T][z-lJ 

A A^ A A^ 

- Izexp(- T)sin(A2- [z-l]}[z-l]‘^[z^I-2zexp(- T)cos(A2- 

+ exp(-A^T)] 

which is 

(A. 1.6) Z{s"^[I 8^4 -Aj^ 8+A2]} - {(I~exp(- T)(cos(A2- 

A *2 .2 

1 - 1/2 1/2 

■ 2 ' 8in(A2- ^)^'^T])z+exp(- -f T)(exp(- T) 

- [co 8(A2~ (A2- -^)"^^^8ln(A2- )z}{ (z-1)"^ 

A A^ 

• [z^I-2zexp(- T)cos(A 2- -^)^'^^T+exp(-A^T) 

where 
(A. 1.7) 

(A. 1.8) A2 = 

Defining 

(A. 1.9) 3 ^ - - ^ _ i 





(A.l. Z{s"^[Is2+Aj^ 8+A2]"^} - {(I-exp(ej^T)tcos(S2T)-3j^B2^8ln(62T))z^ 

+ exp(3j^T(exp(6jT)-[co8(62T)-6j^B2^sin(32T) ])z}{ (z-1) 

' [z^I-2zexp(3j^T)cos(32T)+exp(23j^T) ] } ^ 


I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 

I 


Returning to (A. 1.2), the discrete transfer function T(z) Is then 

(A. 1.11) T(z) - (I-exp (3j^T) [cos (32T)-3j^32^sin ( 32 T) ) z 

+ exp(3j^T) (exp(3j^T)-[cos(32T)-3j^32^sin(32T)) }{z^I-2zexp(3j^T) 
• cos(32T)+exp(23j^T)}"V'^^G 

where 3^^ and 32 defined as In (A. 1.9) and (A. 1.10). 
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ROOTS OP THE SHIFTED LEGENDRE POLYNOMIALS 
AND CHRISTOPPEL WEIGHTS 


1 

1 

1 

I 

1 

I 

I 


ROOTS 

N = 3 

WEIGHTS 

1 .1270166E-1 
5.0000000E-1 
8 . 8729833 E-I 

N * 5 

2.7777777E-1 

i)..4444444E-l 

2.7777777E-1 

4.6910077E-2 

2.3076534E-1 

5.0000000E-1 

7 . 6923465 E-I 

9 . 5308992 E-I 

N = 7 

1.1846344E-1 

2.3931433E-1 

2.8444444E-1 

2.3931433E-1 

1.1846344E-1 

2. 544604 3E-2 

1 . 2923440 E-I 
2.9707742E-1 
5.0000000E>1 
7.0292257E-1 
8 . 7076559 E-I 
9 . 7455595 E-I 

N = 9 

6.4742483E-2 
1 . 3985269E-1 
I. 9091502 E-I 
2.0897959E-1 
I. 9091502 E-I 
1 . 3985269 E-I 
6.4742483E-2 

1 . 591988 OE -2 
8.1984446E-2 
1.9351428E-1 
3 . 5787328 E-I 

5.0000000E-1 

6 . 621267 IE-I 

8.0668571E-1 

9 .I 801555 E-I 

9.8408011E-1 


4.0637194E-2 
9 .O 32408 OE -2 
1 . 3030534E-1 
1.5617353E-1 
I. 6511967 E-I 
1.5617353E-1 
1 . 3030534E-1 
9 .O 32408 OE -2 
4 . 0637194 E -2 


I 
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ROOTS 


WEIGHTS 


N » 11 


1 .O88567OE-2 
5.6468700E-2 
l.?^92399E-l 
2.itOik51Q5E-l 
3.6522842E-1 
5.0000000E-1 
6.3477157E-1 
7.595^806e-1 
8.650760OE-I 
9.^353129E-1 
9.8911432E-I 


N “ 13 


7.9084726E-3 

4.120080OE-2 

9.92l095^E-2 
1 . 7882533E-1 
2.7575362E-1 
3.8477084E-1 
5.0000000E-1 
6.1 52291 5E-1 
7.2424637E-1 
8.2117466E-1 
9.0078904E-1 
9.5879919E-1 
9.9209152E-I 

1 = 15 

6.0037409E-3 
3.1363303E-2 
7.5896708E-2 
1. 377911 3E-1 
2.145139IE-I 
3.0292432E-1 
3. 994029 5E-1 
5.0000000E-1 
6.0059704E-1 
6.9707567E-1 
7.8548608E-1 
8.6220886E-I 
9.2410329E-1 
9.6863669E-I 
9.9399625E-1 


2.7834283E-2 
6.2790184e-2 
9.3145105E-2 
I.I659688E-I 
1.3140227E~1 
1.3646254E-1 
1 .3140227E-1 
I.I659688E-I 
9.3145105E-2 
6.2790184E-2 
2.7834283E-2 


2.0242002E-2 

4.6060749E-2 
6.9436755E-2 
8.907299OE-2 
1.0390802E-1 
1.1 3141 59E-1 
1.1627577E-1 
1.1 3141 59E-1 
I.0390802E-I 
8.907299OE-2 
6.9436755E-2 
4.6060749E-2 

2.0242002E-2 


I.537662OE-2 
3. 516302 3E-2 
5. 357961 OE-2 

6.9785338E-2 
8.3134602E-2 
9.30805OOE-2 
9.9215742E-2 
1.01 2891 2E-1 
9.9215742E-2 
9.30805OOE-2 
8.3134602E-2 
6.9785338E-2 

5. 357961 OE-2 

3.5183023E-2 

I.53766OOE-2 



I 

1 

j NEGATIVES OP LOGARITHMS OP ROOTS OP SHIPTED 

LEGENDRE POLYNOMIALS 



N = 3 

N = 11 

N - 15 

•* 

2.183011 

4.520308 

5 . 11537 i' 


0 . 5931^7 

0.119574 

2.874069 

2.003044 

3 . 46211 r 
2.578382 

- 

1.425235 

1.982016 


N = 5 

1.007232 

0.693147 

0.454490 

1 . 539381 

1.194272 


3.059523 

0.917764 


1 . 46635 ^ 

0.693147 

0.275032 

0.693147 


0.144938 

0.509831 


0.262359 

0 .o 48 o 46 . 

0.058126 

0.360861 


0.010945 

0.241453 


N = 7 

0.148258 

0.078931 

N = 13 



0.031866 


3.671195 

4.839821 

0.006022 


2.046127 

3.189298 



1.213762 

2.310507 



0.693147 

1.721346 



0.352509 

1 . 288247 



0.138382 

0.955107 


T 

0.025775 

0.693147 



N = 9 

0.485760 

0.322624 


] 

4.140187 

0.197019 

0.104484 



2.501226 

0.042074 


i 

1.643438 

1.085084 

0.007940 


0.693147 

0.412298 

0.214821 

O.O05541 

0.016048 




I 


rj r 


I 

I 

I 

1 




irrwt. r’-Yrga»'<g«ara>--g--4->*a>^ 


Appradlx 3 AX'b Version of Golub Algorithm 


SUbROUTINE CJLUB(A.Xflt.MfN) 

C A(hfH) f 9(H) GIVEN WITH M)N SOLVES FOR' X< N ) SnrM tha' 

C J IS A HIHinUH 

C HETHOD or C. GOLUP* NUHERISCHE HATHEHATIK 7. 206- 
IMPLICIT DOUIH.E PRECISION (ID 
REAL A(10>10)»X(10)»B(10)»U(50) 

C PERFORH N ORTHOGONAL TRANSFORMATIONS TO A^ . , . ) TO 

C UPPER TRIANGULAPIZE THE MATRIX 

DO 3010 Kri,N 
DSUM-O.ODO 
DO 1010 I«KrM 
DAJ«A( I»K ) 

1010 DSUMaDSUM^DAJ»«2 
DAI*A(KfK ) 

DSIGMA=DSICN( DSQRT( DSUM )fDAI ) 

DBI=DSORT( l.ODO+DAI/DSIGMA) 

U( K >=DBI 
FACT»DFACT 
KPLU5*K+1 
DO 1020 I«KPLUSfM 
1020 U( I )«FACT*A( IfK) 

C I-U*U»*T IS SYHMETRlCf ORTHOGONAL MATRIX WHICH WHCM 

C TO A(.».) WILL ANNIHILATE THE ELEMENTS BELOW THE DlAC-nwAL 

DO 2030 J=K»N 

C APPLY THE ORTHOGONAL TRANSFORMATION 

r hL’T= 0.0 
DO 2010 I*K»M 
2010 FACT=FACT+U( I >»A( I f J ) 

DO 2020 I=K»M 

2020 A< I » J )»A( I » J )-FACTiU( I ) 

2030 CONTINUE 
FACT»0.0 
DO 2040 I»K»M 
2040 FACT=FACT-»-U( I >*B( I ) 

DO 2050 I*K.M 
2CC0 I )»B( I )-FACT*U< I ) 

3010 CONTINUE 

C BACK SUBSTITUTE TO RECURSIVELY YIELD X( . > 

X(N)BB(N)/A(NrN) 

LIM=N--1 

DO 4020 I«lfLIH 
IROU»N-I 
SUH=0.0 
DO 4010 J«..i 

4010 SUM=SUM+X( N-J+1 >*A( IROUrN-J+1 ) 

4020 X( IROW )»( B< IROW >-SUM )/A( IROWr IROW > 

RETURN 

END 


From Clsrebout , [ 30 ] . 
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